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ABSTRACT 

Pulsar timing arrays (PTAs) are designed to detect gravitational waves with periods 
from several months to several years, e.g. those produced by by wide supermassive 
black- hole binaries in the centers of distant galaxies. Here we show that PTAs are 
also sensitive to mergers of supermassive black holes. While these mergers occur on 
a timescale too short to be resolvable by a PTA, they generate a change of metric 
due to non-linear gravitational-wave memory which persists for the duration of the 
experiment and could be detected. We develop the theory of the single-source detec- 
tion by PTAs, and derive the sensitivity of PTAs to the gravitational-wave memory 
jumps. We show that mergers of lO^M© black holes are 2 — a-detectable (in a direc- 
tion, polarization, and time-dependent way) out to co-moving distances of ~ 1 billion 
light years. Modern prediction for black-hole merger rates imply marginal to modest 
chance of an individual jump detection by currently developed PTAs. The sensitivity 
is expected to be somewhat higher for futuristic PTA experiments with SKA. 



1 INTRODUCTION 



Bursts of gravitational waves leave a permanent imprint 
on spacetime by causing a small permanent change of 
the me t ric, as computed in the tran s verse traceless gauge 
f Payne" '1983'; ' Christodouloul Il99ll : iBlanchet fc DamouJ 
11992 ; _Thornc 1992). This gravitational- wave "memory 
jumps" are particularly significant in the case of merger 
of a b i nary black hole, as was recently pointed out by 
iFavatal (|2009l . hereafter F09). Favata has shown (see Fig- 
ure 1 of F09) that for the case of an equal-mass binary, 
a metric memory jump Sh was of the order of ~ 5 per- 
cent of M/R, where M is the mass of the binary compo- 
nent and R is the co-moving distance to the binary mea- 
sured at redshift (hereafter M is expressed in the geo- 
metric units, i.e. M — GM/(?). Furthermore, Favata has 
argued that the memory jumps were potentially detectable 
by LISA with high signal-to-noise ratio. Favata's memory 
calculations make use of an approximate analytical treat- 
ment of the mergers, and need to be followed up with more 
definitive numerical calculations. Nevertheless, a number of 
analytical models explored in F09 show that the effect is 
clearly of high importance, and thus further investigations 
of detectability of the memory jumps are warranted. 

Recently, there has been a renewed effort to mea- 
sure gravitational waves from widely separated supermas- 
sive black-hole (SMBH) bina ries by using prec i se timing of 
galac tic millisecond pulsars (jjenet et al.ll2005l ; iManchesteJ 
I2OO6I ). In this paper we investigate whether pulsar timing 
arrays (PTAs) could be sensitive to the memory jumps from 
physical mergers of the SMBHs at the end of the binary 's 
life. We demonstrate that modern PTAs l|Manchestedl2006l ). 
after 10 years of operation, will be sensitive to mergers of 
10^ Mq black holes out to ~billion light years; however the 



chances of actual detection are small. Futuristic PTA exper- 
iments, like those pe rformed on the Square Kilometer Array 
dCordes et al.|[2005l l. offer a somewhat better prospect for 
the direct detection of gravitational-wave memory jumps. 



2 THE SIGNAL 

The gravitational waveform from a merger of SMBH pair 
consists of an ac-part and a rfc-part; see Figure 1 of F09. 
The ac-part is short-period and short-lived, and hence is 
undetectable by a PTA. The dc-part is the gravitational- 
wave memory; it grows rapidly during the merger, on the 



timescale of ~ 10M(1 -I- z) 



1O*(M/1O*'M0)(1 + z) 



where M is the mass of the SMBHs (assumed equal) and 
z is the redshift of the merger. After the burst passes, the 
change in metric persists, and as we explain below, it is 
this durable change in the metric that makes the main im- 
pact on the timing residuals. Realistic PTA programs are 
designed to clock e ach of the pulsars with ~ 2-week intervals 
(|Manchesteill200^ . Bailes, private communications). There- 
fore, for M = 10* M0 SMBHs the growth of the memory- 
related metric change is not time-resolved by the timing 
measurements. Moreover, even for M = 10^°A/o SMBHs 
this growth occurs on the timescale much shorter than the 
duration on the experiment. We are therefore warranted to 
treat the rfc-part of the gravitational wave as a discontinuous 
jump propagating through space, 

h{r,t)^haxQ[{t-to)-n-^, (1) 

where ho is the amplitude of the jump, of the order of 
0.05M/i?, 0(t) is the Heavyside function, to is the moment 
of time when the gravitational-wave burst passes an ob- 
server, r is the location in space relative to the observer. 
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and n is the unit vector pointed in the direction of the wave 
propagation. Here and below we set c = 1. We have used 
the plane- wave approximation, which is justified for treating 
extragalactic gravitational waves as they propagate through 
the Galaxy. 

For a single pulsar, the frequency of the pulse- 
arrival I' responds to a pla ne gravitational wave accord- 
ing to the following equation l|Estabrook fc Wahlguistll 19751 : 
iHellings fc Downs, 198^ : 



V 

where 



B{e, (j)) X [hit) --h(t~r-r cos 9)] , 



-cos(2<^)(l-cos6'). 



(2) 



(3) 



Here r is the Earth-pulsar distance at an angle 6 to the 
direction of the wave propagation, (fi is the angle between the 
wave's principle polarization and the projection of the pulsar 
onto the plane perpendicular to the propagation direction, 
and h{t) is the gravitational-wave strain at the observer's 
location. Substituting Eq. ([1]) into the above equation, we 
obtain the mathematical form of the signal: 



5u{t) 



= hoB{9, 0) X [e{t - to) - Q{t - ti)] , 



(4) 



where = to + ''(1 + cos 6*). Thus the memory jump would 
cause a pair of pulse frequency jumps of equal magnitude 
and the opposite sign, separated by the time interval r(l -|- 
cosO). Since typical PTA pulsars are at least ~ 10^ light 
years away, a single merger could generate at most one of 
the frequency jumps as seen during the ~ 10 years of a 
PTA experiment. The timing residuals from a single jump 
at t = to are given by 



m{t) = B{e, (p)ho X e(t - to) X (t - to). 



(5) 



For a single pulsar the frequency jump is indistinguish- 
able from a fast glitch, and therefore single-pulsar data can 
only be used for placing upper limits on gravitational-wave 
memory jumps. The situation would be different for an ar- 
ray of pulsars, where simultaneous pulse frequency jumps 
would occur in all of them at the time t = to when the 
gravitational-wave burst would reach the Earth. Therefore 
a PTA could in principle be used to to detect memory jumps. 



3 SINGLE-SOURCE DETECTION BY PTAS. 

In this section we develop a mathematical framework 
for the single-source detection by a PTA. Our formalism 
is essentially Bayes i an an d follows closely the spirit of 
Ivan Haasteren et al.l l|2009l . hereafter vHLML), although we 
will make a connection with the frequentist Wiener-filter es- 
timator. We will then apply our general formalism to the 
memory jumps. The reader uninterested in mathematical 
details should skip the following subsection and go straight 
to the results in section O 

There is a large body of literature on the single-source 
detection in the gravitational- wave community (.Finn, ,1993 : 
IOwenlll996l : iBradv et al.lll998l) . The techniques which have 
been developed so far are designed specifically for the inter- 
ferometric gravitational- wave detectors like LIGO and LISA. 
There are several important modifications which need to be 



considered when applying these techniques to PTAs, among 
them 

1. Discreteness of the data set. A single timing residual 
per observed pulsar is obtained during the observing run; 
these runs are separated by at least several weeks. This is in 
contrast to the continuous (for all practical purposes) data 
stream in LIGO and LISA. 

2. Subtraction of the systematic corrections. The most 
essential of these is the quadratic component of the timing 
residuals due to pulsar spindown, but there may be others, 
e.g. jumps of the zero point due to equipment change, an- 
nual modulations, etc. 

3. Duration of the signal may be comparable to the 
duration of the experiment. This is the case for both 
cosmological stochastic background considered in vHLML, 
and for the memory jumps considered here. Thus frequency 
domain methods are not optimal, and time-domain formal- 
ism should be developed instead. 

The Bayesian time-domain approach developed in 
vHLML and in this subsection is designed to tackle these 
3 complications. 

Consider a collection of A'^ timing residuals Stp obtained 
from clocking a number of pulsars. Here p is the composite 
index meant to indicate both the pulsar and the observing 
run together. Mathematically, we represent the residuals as 
follows: 



5tp = Ax s{tp) + Stp + Q{tp). 



(6) 



Here s{tp) and A are the known functional form and un- 
known amplitude of a gravitational- wave signal from a single 
source, Stp is the stochastic contribution from a combination 
of the timing and receiver noises, and 

Q{ip) ~ ^m,^mfm{tp) (7) 

is the contribution from systematic errors of known func- 
tional forms fmitp) but a- priory unknown magnitudes ^m. 
Below we shall specify Q{tp) to be the unsubtracted part 
of the quadratic spindown, however for now we prefer to 
keep the discussion as general as possible. We follow van 
Haasteren et al. (2009) and rewrite Eq. (0) in a vector form: 



St = As + St + F^. 



(8) 



Here the components of the column vectors St, St , s, and ^ 
are given by Stp, Stp, s{tp), and ^,n, and _F is a non-square 
matrix with the elements Fpm = fmitp). Henceforth we as- 
sume that 5tp is the random Gaussian process, with the 
symmetric positive-definite coherence matrix C: 

Cp, =< St';St':^ > . (9) 

We can now write down the joint probability distribution 
for A and 



P{A,£,m\St) 



{l/M)Po(A,£,m) X 



(10) 



exp 



-]^(5t- As - Fi)'^ X C"^x, 



[Si- As- F^)] 



Here Po{A,^rn) is the prior probability distribution, and M 
is the overall normalization factor. We now assume a fiat 
prior Po{A, Lm) — const, and marginalize over ^ in precisely 
the same way as shown in the Appendix of vHLML. As a 
result, we get the following Gaussian probability distribution 
for A: 
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P{A\St) = 



27ro- 



exp 



2cr2 



(11) 



Here, the mean value A and the standard deviation a are 
given by 



A 



and 

where 
C' = C" 



1/2 



(12) 

(13) 
(14) 



It is instructive and useful to re-write the above equations 
by introducing an inner product (x, y) defined as 

{x,y) = x^C~^y. (15) 

Let us choose an orthonormal basi£] fi in the subspace 
spanned by /„, so that {fijj) 
a projection operator 



We also introduce 



R 



fm){frr 



(16) 



so that Rx — X — Timifm, x) fm- AH the usual identities 
for projection operators are satisfied, i.e. R^ = R and 
{Rx, Ry) = (x, Ry) . We can then write 



A 



,RSt) 



{s,Rs) ' 
and 

a = {s,Rs)-^^\ 



(17) 



(18) 



If there are no systematic errors that need to be removed, 
than R = 1 and the Eqs ((17} and (|18|) represent the time- 
domain version of the Wiener-filter estimator. 



3.1 Other parameters 

So far we have assumed that the gravitational-wave signal 
has a known functional form but unknown amplitude, and 
have explained how to measure or constrain this amplitude. 
In reality, the waveform s{ff,tp) will depend on a number of 
a-priori unknown parameters ff, such as the starting time of 
the gravitational-wave burst and the direction from which 
this burst has come. These parameters enter into the proba- 
bility distribution function through s in Eq. (|10p , and gener- 
ally their distribution functions have to be estimated numer- 
ically. The e stimates can be done via the matched filtering 
l|Owenlll996l ) or by performing Markov-Chain Monte-Carlo 
(MCMC) simulations. In section [S] we will demonstrate an 
example of an MCMC simulation for the memory jump. In 
this section, we show how to estimate an average statistical 
error on ff for signals with high signal-to-noise ratios. 

Let us begin with a joint likelihood function for the 
amplitude A and other parameters 77: 

L(A, 77) = -{l/2){As{ff) - St, R{As{ff) - St)) + Const. (19) 



We now fix A to its maximum-likelihood value 
{s{ff) , RSt) / {s{ff) , Rs{ff)) , and average over a large number 
of statistical realizations of the noise St . The so-averaged 
likelihood function is given by 

(1/2)A? 



iav (ff) = 



{s{fi),Rm} 

{s{fft),Rs{ff))^ 



Mfft),Rs{fft)){siff),Rs{ff))- 



(20) 



where At,fjt are the true values for the signal present in all 
data realizations. We have omitted the additive constant. 

The expression in the square bracket is positive-definite, 
and Lav is quadratic in rf— fft for the values of ff close to the 
true values, 



Lav(r;)~-(l/2)(rf-rfOG(»7-rfO, 



(21) 



where G is the positive-definite Fisher information matrix. 
Its elements can be expressed as 



A!/{s,Rs) 



drji drjj 



{S,R^){S,R^) 
dm drjj 



(22) 



evaluated at rj — -qt. The inverse of G specifies the average 
error with which parameters ff can be estimated from the 
data. 



4 DETECTABILITY OF MEMORY JUMPS 

We now make an analytical estimate for detectability of the 
memory jumps. For simplicity, we assume that all of the pul- 
sar observations are performed regularly so that the timing- 
residual measurements are separated by a fixed time At, 
and that the whole experiment lasts over the time inter- 
val [~T, T] (expressed in this way for mathematical conve- 
nience). Furthermore, we assume that the timing/receiver 
noise is white, i.e. that for a pulsar a 



{St" St]) = alSt. 



(23) 



This assumption is probably not valid for some of the mil- 
lisecond pulsars (Verbiest et al., in prep., van Haasteren et 
al., in prep.). We postpone discussion of the non-white noises 
to future work. 

To keep our exposition transparent, we consider the case 
when the array consists of a single pulsar a; generalization 
to several pulsars is straightforward and is shown later this 
section. Finally, we assume that the systematic error Q{ti) 
comprises only an unsubtracted component of the quadratic 
spindown. 



Q{t,) = Ao + Alt, + Aitt 



(24) 



where Aq, Ai, and A2 are a-priori unknown parameters. 

We now come back to the formalism developed in the 
previous section. The inner product defined in Eq. (|15|l takes 
a simple form: 



,y) = -^y2x{ti)y{ti) 



(25) 



^ This is always possible by e.g. the Gramm-Schmidt procedure. 



x[t)y(t)dt. 
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where we have assumed At <^ T and have substituted the 
sum with the integral in the last equation. We now choose or- 
thonormal basis vectors /i,2,3(t) which span the linear space 
of quadratic functions: 



(26) 




From Eq. ([5]) the gravitational-wave induced timing residu- 
als are given by 5t{t) — h()s{t), where 



s{t) = B{e, (p) X e{t - to) X (t - to)- 



(27) 



The expected measurement error of the jump amplitude ho 
is given by Eq. (fT8|l : 



1=1,2,3 



-1/2 



{s,s) 



Substituting Eqs. ((23, and ^ into Eq. 
after some algebra: 



(28) 
one gets 



B{e, </>) T y iVp3 (i - 



(29) 



Here = 2T/At is the number of measurements, and 

p=l- [to/Tf . (30) 

For an array consisting of multiple pulsars, and with the 
assumption that the timing residuals are obtained for all 
of them during each of the A'^ observing runs, the above 
expression for is modified as follows: 



T yivp3(i-i|p) 



(31) 



same amplitude of timing/receiver noise Ua ~ cr, the a^g in 
Eq. (I3ip is given by 



(34) 



4. While the timing precision of future timing arrays is 
somewhat uncertain, it is instructive to consider a numerical 
example. Lets assume T = Syr (i.e., the 10-year duration of 
the experiment), = 250 (i.e., roughly bi-weekly timing- 
residual measurements), A'j, = 20 isotropically-distributed 
pulsars (this is the current number of clocked millisec- 
ond pulsars), and aa ~ 100ns (this sensitivity is currently 
achieved for only several pulsars). Then for optimal arrival 
time to = ±r/\/5, the array sensitivity is 

(Tho = 4.5 X 10"^^ (35) 

For a binary consisting of two black holes of the mass M, 
the memory jump is estimated in F09 to be 



' R 0.05 \10»Mq 



10^ light-years 
R 



,(36) 



where -q ~ 0.05 is the direction-dependent numerical param- 
eter. In this example, the pulsar-timing array is sensitive 
to the memory jumps from black-hole mergers at redshifts 
z < zo, where zo ~ 0.1 for M = IO^A/q, and zo ~ 1 for 
M = IO'^Mq. 



4.1 Arrival time 

It is possible to estimate the array's sensitivity to the 
memory-jump arrival time, to. We use Eqs. ((22} and (|27|) . 
and after some algebr£0 get 



o-to ■- 

where 
1 

x^(p) 



(38) 



f 1 + - 2p) - 3(l-p)p [l-(5/8)p]-^ 
2^ V r V 2[1 - (15/16)p] ^ ' 



where 



Cefi = 



-1/2 



(32) 



Several remarks are in order: 

1. The error aha diverges when p — 0, i.e. when to = 
±T. This is as expected: when the memory jump arrives at 
the beginning or at the end of the timing-array experiment, 
it gets entirely fitted out when the pulsar spin frequency is 
determined, and is thus undetectable. 

2. Naively, one may expect the optimal sensitivity when 
the jump arrives exactly in the middle of the experiment's 
time interval, i.e. when to — 0. This is not so; the optimal 
sensitivity is achieved for to/T — ±l/y/E when the error 
equals 



(Jho 



(33) 



3. The sky-average value for B'\e,(j>) is 1/6. Therefore, 
for an array consisting of a large number of pulsars Np which 
are distributed in the sky isotropically and which have the 



4.2 Source position 

The array's sensitivity gravitational-wave memory is depen- 
dent on source position since the number and the position 
of the pulsars in current PTAs is not sufficient to justify 
the assumption of isotropy made in Eq. (I34|l . We will there- 
fore calculate the value of [X^^ {Qa,4'a)\ ^^'^ for current 
PTAs. Since the polarisation of the gravitational-wave mem- 
ory signal is an unknown independent parameter, we average 
over the polarisation and obtain for the angular sensitivity: 



-1/2 



^ A useful identity: 

a[(t-to)e(f-fo)] 

dto 



^i(l+cos^ 



-e(t-to). 



-1/2 



(40) 



(41) 



(37) 
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Figure 1. The relative sensitivity cr^ff for the pulsars of the Eu- 
ropean Pulsar Timing Array. The scaling has been chosen such 
that a value of 1 indicates that the same sensitivity for that source 
position would have been achieved with a perfect isotropic PTA 
(i.e. B2 = 1) 




Figure 2. The relative sensitivity a^g for the pulsars of the 
Parkes Pulsar Timing Array. The scaling has been chosen such 
that a value of 1 indicates that the same sensitivity for that source 
position would have been achieved with a perfect isotropic PTA 
(i.e. S2 = 1) 

Here we have assumed that all pulsars have equal timing pre- 
cision. 03 and 8s are the position angles of the gravitational- 
wave memory source, and 8a is the polar angle of pulsar a 
in a coordinate system with {(l)a,6a) at the north-pole. In 
figure [1] and [2] the sensitivity to different gravitational- wave 
memory source positions is shown for respectively the Eu- 
ropean Pulsar Timing Array and the Parkes Pulsar Timing 
Array projects. 



5 TESTS USING MOCK DATA 

We test the array's sensitivity to gravitational-wave mem- 
ory signals using mock timing residuals for a number of mil- 
lisecond pulsars. In this whole section, all the mock timing 



residuals were generated in two steps: 

1) A set of timing resid uals was generated using the pulsar 
timing package tempo2 (|Hobbs et al.ll200^ ). We assume that 
the observations are taken tri-weekly over a time-span of 10 
years. The pulsar timing noise was set to 100 ns white noise. 

2) A gravitational-wave memory signal was added according 
to Eq. ((5]), with a memory-jump arrival time set to be opti- 
mal for sensitivity: — The direction and polarisation 
of the gravitational-wave memory signal were chosen ran- 
domly - the coordinates happened to have declination 90°. 
In the following subsections we describe tests which have 
fixed parameters for step 1, but systematically varied am- 
plitude for step 2, and we use these tests to study the sen- 
sitivity of the array. 

5.1 Used models 

In principle, we would like to realistically extrapolate the re- 
sults we obtain here for mock datasets to future real datasets 
from PTA projects. Several practical notes are in order to 
justify the models we use here to analyse the mock datasets: 

1) From equation (|24|) onward, we assume that the 
systematic-error contributions to the timing residuals con- 
sist only of the quadratic spindown. In reality, pulsar ob- 
servers must fit many model parameters to the data, and 
have developed appropriate fitting routines within timing 
packages like tempo2. Similar to the quadratic spindown dis- 
cussed in this paper, all the parameters of the timing model 
are linear or linearised in tempo2, and therefore those pa- 
rameters are of known functional form. Since the subtraction 
of quadratic spindown decreases the sensitivity of the PTA 
to gravitational-wave memory signals, we would expect the 
same thing to be true for the rest of the timing model. 

2) The error-bars on the pulse arrival time obtained from 
correlating the measured pulse-profile with the template of 
the pulse-profile are generally not completely trusted. Many 
pulsar astronomers invoke an extra "fudge" factor that ad- 
justs the error-bars on the timing-residuals to make sure that 
the errors one gets on the parameters of the timing-model 
are not underestimated. Usually the "fudge" factor, which 
is known as an efac value is set to the value which makes 
the reduced of the timing solution to be equal to 1. 

In order to check the significance of both limitations 1 
& 2, we perform the following test. We take a realistic set 
of pulsars with realistic timing models: the pulsar positions 
and timing models of the PPTA pulsars. We then simulate 
white timing-residuals and a gravitational-wave memory sig- 
nal with amplitude ho = 10~^^, and we produce the poste- 
rior distribution of Eq. (|ll|l in three different ways: 

a) We marginalise over only the quadratic functions of 
Eq. ([26} , which should yield the result of Eq. (|3TJ . 

b) We marginalise over the all timing model parameters in- 
cluded in the tempo2 analysis when producing the timing- 
residuals. 

c) We marginalise over all the timing model parameters, 
and we also marginalise over the efac values using the nu- 
merical techniques of vHLML. By estimating the efac value 
simultaneously with the gravitational-wave memory signal, 
we are able to completely separate the two effects. Note that 
this procedure will not destroy information about the rela- 
tive size of the error-bars for timing-residuals of the same 
pulsar. 
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Sensitivity to GW-memory 



Marginalise over quadratics 
Ivlarginaiise over compiete timing-modei 
lylarginaiise over timing model and etacs 




-4e-15 -3e-15 -2e-15 -1e-15 1e-15 2e-15 3e-15 4e-15 5e 

Amplitude ot GW-memory signai 

Figure 3. The posterior distribution of tlie gravitational-wave 
memory amplitude. Tlie two solid lines are the result of an anal- 
ysis where we only analytically marginalise over the full timing 
model or just quadratic spindown. The points are the result of a 
marginalisation over the full timing model and the efac values as 
well. From the Gaussians, the sensitivity can be reliably estimated 



lylarginaiised posterior distribution tor GW-memory signal 




Starting time signal [yr] 



Figure 4. The marginalised posterior distribution for the 
gravitational-wave memory signal amplitude and arrival time of 
the jump. In this case a dataset was analysed that did not contain 
any gravitational-wave memory signal. 



at: 



6.5 X 10- 



We present the result of this analysis in Figure [3] Based 
on the 185 observations per pulsar in the dataset and the di- 
rection of the gravitational- wave memory signal, we can cal- 
culate the theoretical sensitivity of the array using Eq. (|33p 
and Eq. p4t . This yields a value of: 

o-ho = 6.4 X 10"'^ (42) 

We can also calculate this value for the three graphs in Fig- 
ure [3] The three graphs lie close enough on top of each other 
to conclude that one value applies to all three of them: 

o-hQ = 6.5 X 10-'^ (43) 

which is in good agreement with the theoretical value. It 
appears that both note 1 and 2 mentioned above are not of 
great influence to the sensitivity of PTAs to gravitational- 
wave memory detection; the theoretical calculations of this 
paper are a good representation of the models mentioned in 
this section. 

5.2 Upper-limits and detecting the signal 

When there is no detectable gravitational- wave memory sig- 
nal present in the data, we can set some upper-limit on the 
signal amplitude using the algorithm presented in this paper. 
Here we will analyse datasets with no or no fully detectable 
gravitational- wave memory signal in it, and a dataset with a 
well-detectable signal using the MCMC method of vHLML. 
We will calculate the marginalised posterior distributions 
for the 5 parameters of the gravitational-wave memory sig- 
nal. The interesting parameters in the case of an upper- 
limit are the amplitude and the arrival time of the jump. 
A marginalised posterior for those two parameters are then 
presented as two-dimensional posterior plots. Note that the 
difference with the analysis in section 15.11 is that we vary 
all gravitational-wave memory parameters, instead of only 
the amplitude. Note that we do marginalise over all the efac 
values as discussed in section [5TT1 unless stated otherwise. 



In Figure|4]we show the result of an analysis of a dataset 
where we have not added any gravitational-wave memory 
signal to the timing-residuals. The 3 — a contour is drawn, 
which serves as an upper-limit to the memory amplitude. 
We see that we can exclude a gravitational-wave memory 
signal at ^ = of amplitude 3 x 10~^^ and higher. We 
see that this value is over a factor of 4 higher than what is 
predicted by Eq. (|42[l . This is to be expected, since: 

1) We give a 3 — (T limit here, instead of the 1 — a sensitivity. 

2) We also marginalise over the arrival time and other pa- 
rameters of the memory signal, reducing the sensitivity. 
Because of these reasons, we argue that the minimal upper- 
limit one can set on the gravitational-wave memory sig- 
nal using a specific PTA is the sensitivity calculated using 
Eq. (|42j muhiplied by 4. 

Next we produce a set of timing-residuals with a mem- 
ory signal of amplitude ho = 10~^^. According the result 
mentioned above, the memory signal should not be resolv- 
able with this timing precision. The result is shown in Fig- 
ure [5l We see that we can indeed merely set an upper- limit 
again. In order to check the effect of marginalising over the 
efac values as mentioned in section ISTTI we also perform an 
analysis where we pretend we do know the efac values prior 
to the analysis. The result is shown in Figure |6l We see no 
significant difference between the two models. 

Finally, we also analyse a dataset with a gravitational- 
wave memory signal with an amplitude larger than the 3 x 
10~^^ upper- limit of the white set mentioned above. Here we 
have added a memory signal with an amplitude of 10~^*. In 
Figure [7] we see that we have a definite detection of the 
signal: if we consider the 3 — ct contours, we see that we can 
restrict the gravitational-wave memory amplitude between 
[6.6 x 10"^'^, 1.35 X 10""] . Again, this value is higher than 
the value predicted by Eq. (|42|l due to us including more 
parameters in the model than just the memory amplitude. In 
Figure |8] we see that we can also reliably resolve the position 
of the source in this case. 
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Marginalised posterior distribution for GW-memory signal 
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Figure 5. The marginalised posterior distribution for the 
gravitational-wave memory signal amplitude and arrival time of 
the jump. Here a gravitational-wave signal with an amplitude of 
10~^^ was added to the white residuals. The contour drawn is 
the 3 — (T contour. 



Figure 7. The marginalised posterior distribution for the 
gravitational-wave memory signal amplitude and arrival time of 
the jump. Here a gravitational-wave signal with an amplitude of 
10"^* was added to the white residuals, indicated with a in 
the figure. The contours drawn are the 1 — ct and 3 — ct contours. 



Ivlarginaiised posterior distribution tor GW-memory signal 





starting time signal [yr] 



Figure 6. The marginalised posterior distribution for the 
gravitational-wave memory signal amplitude and arrival time of 
the jump. Here a gravitational-wave signal with an amplitude of 
10~^^ was added to the white residuals. This analysis has been 
done without marginalising over the efac values. The contour 
drawn is the 3 — cr contour. 



Figure 8. The marginalised posterior distribution for the sky 
location of the gravitational-wave memory signal. We can see here 
that we can marginally determine the direction of the source. The 
source positions used to generate the residuals were: (declination, 
right ascension) = (90", 12.4hr). 



6 DISCUSSION 

In this paper, we have shown that gravitational-wave mem- 
ory signals from SMBH binary mergers are in principle de- 
tectable by PTAs, and that 2 ~ a constraints are possible 
on M = 10* Mq mergers out to redshift of ~ 0.1 (while 
those with M = 10^" Mq should be detectable throughout 
the Universe) . How frequently do these mergers occur during 
the PTA lifetime? Recent calculations of lSesana et al.1 l|2007l . 
SVH) are not too encouraging. SVH compute, for several 
models of SMBH merger trees, the rate of SMBH mergers 
as seen on Earth, as a function of mass (their figure Id), as 
well as a multidute of other parameters for these mergers. 
From their plots one infers few x 10"^ -10"^ PTA-observable 
mergers per year, which converts to at most 0.1 — 0.01 de- 
tected mergers during the PTA lifetime of ~ 10 years (NB: 



during the PTA existence, only a fraction of time will spent 
near the arrival times with optimal sensitivity). It is conceiv- 
able that SVH estimates are on the conservative side, since 
the mergers of heavy black holes may be stalled (due to the 
"last parsec" problem) and may occur at a significantly later 
time than the mergers of their host halos. In this case, some 
fraction of high-redshift mergers may be pushed towards 
lower redshifts and become PTA-detectable. Detailed cal- 
culations are needed to find out whether this process could 
substantially increase the rate of PTA-detectable mergers. 
It is also worth pointing out that a futuristic PTA experi- 
ment based on a Square Kilometer Array may attain up to 
an order of magnitude higher sensitivity that the currently 
developed PTAs. 

The methods presented in this paper are useful be- 
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yond the particular application that we discuss. The algo- 
rithm presented here is suitable for any single-source detec- 
tion in general when the gravitational waveform has known 
functional form. Further applications will be presented else- 
where. 

6.1 comparison with other work 

When this paper w as already finished, a preprint by 
IPshirkov erall l|2009l , PBP) has appeared on the arxiv which 
has carried out a similar analysis to the one presented here. 
Our expressions for the signal-to-noise ratio for the mem- 
ory jump agree for the case of the white pulsar noise. PBPs 
treatment of cosmology is more detailed than ours, while the 
moderately pessimistic predicted detection rates are broadly 
consistent between the 2 papers. Our method for signal ex- 
traction is more generally applicable than PBS's since it is 
optimized for any spectral type of pulsar noise, takes into 
consideration not just the signal magnitude but also other 
signal parameters, and is tested on mock data. 
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